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Lawrence  Radiation  Laboratory,  University  of  California 
Livermore,  California 

April  19,  1963 

INTRODUCTION 

The  first  requirement  in  the  calculation  of  elastic -plastic  flow  is  a 
formulation  of  the  equation  of  state.  The  equation  of  state  must  describe 
elastic,  elastic-plastic,  and  hydrodynamic  flow.  The  appropriate  yield 
criteria  must  be  included  in  the  latter  two  regimes.  The  literature  includes 
many  complicated  forms  of  equations  of  state,  some  of  which  have  been 
developed  to  aid  the  mathematics  in  the  analytic  solution  of  the  equations  of 
motion.  However,  since  numerical  techniques  will  be  considered  here,  the 
equations  of  motion  are  completely  independent  of  any  rheological  equation  of 
state  and  any  form  may  be  used.  The  object  of  the  equation  of  state  will  be 
to  provide  a  theoretical  description  applicable  to  a  wide  class  of  practical 
problems,  but  using  simple  idealizations  of  the  outstanding  features  of  the 
real  phenomenon. 

The  problems  of  greatest  present  interest  pertain  to  metal  plasticity, 
but  detailed  description  of  rate- dependent  processes,  for  example,  are  still 
not  well  enough  defined  experimentally.  Therefore,  the  plastic  state  will  be 
described  by  continuously  adjusting  the  stresses  such  that  the  yield  strength 
of  the  material  is  not  exceeded.  More  sophisticated  descriptions  may  be 
included  as  they  seem  indicated  by  experiment. 

This  article  is  arranged  in  three  parts”  Part  I  -  Equation  of  State, 

Part  II  -  One-Dimensional  Elastic -Plastic  Flow,  and  Part  III  -  Two- 
Dimensional  Elastic-Plastic  Flow. 
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PART  I.  EQUATIONS  OF  STATE 


A.  Elastic  Region 


We  are  only  considering  media  which  have  the  same  material  properties 

in  all  directions  (isotropic  media). 

In  x,  y,  z  coordinates  the  state  of  stress  in  a  continuous  media  is 

defined  at  a  given  point  by  six  stress  components,  o^.,  cry,  o^,  t^,  t^, 

and  t  (Ref.  1,  p.  14)-  It  is  always  possible  to  choose  coordinate  axes  such 

xy 

that  the  shear  stresses  at  a  given  point  are  zero  i.  e.  ,  such  that  t  ~  Tzx 

=  t  =0  (Ref.  2,  p.  215).  Any  three  orthogonal  axes  such  that  the  above 
xy 

condition  results  are  called  principal  axes  for  the  point  considered.  The 
stresses  in  the  directions  of  the  principal  axes  on  surfaces  normal  to  these 
axes  are  called  principal  stresses,  denoted  by  or  j  ,  or^ ,  and  cr^. 

A  perfectly  elastic  material  is  characterized  by  a  linear  correspondence 
between  stress  and  strain.  Hooke's  law  is  used  to  describe  the  stress  at  a 
point  resulting  from  a  strain  at  this  point.  The  strain  itself  results  from  a 
force  displacing  particles  in  the  media.  Hooke's  law  in  terms  of  an 
incremental  strain  resulting  in  an  incremental  stress  may  be  written  as. 

V  ,  _  . 

o-j  =  y  +  2pe j  , 

<r 2  -  y  +  2^2’  ^ 


cr3=  y+2pi3. 


Here  X  and  jjl  are  the  Lame  constants  ,  and  and  are  the  strain 

rates  in  the  direction  given  by  the  subscripts;  V  =  volume. 

The  dot  means  a  time  derivative  along  a  particle  path.  It  must  be  noted 
that  the  time  derivative  provides  a  desired  ordered  sequence  for  the  incre¬ 
mental  stress -strain  relationship,  but  this  does  not  mean  that  a  rate- 

3 

dependent  st res s  -  strain  relationship  has  been  introduced.  Hooke's  law  used 
in  this  way  gives  natural  strain,  which  means  that  the  strain  of  an  element  is 
referred  to  the  current  configuration  instead  of  the  original  configuration. 

The  stress  behavior  of  a  material  can  be  thought  of  as  being  composed 
of  a  stress  associated  with  a  uniform  hydrostatic  pressure  (all  three  normal 
stresses  equal)  plus  a  stress  associated  with  the  resistance  of  the  material 
to  shear  distortion.  In  describing  yielding  and  plastic  flow,  we  will  want  to 
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limit  only  the  stress  contributions  that  are  due  to  shear  distortions.  There¬ 
fore,  we  will  decompose  each  of  the  stresses  (Tj  ,  cr^,  and  o*3  into  a  hydro¬ 
static  component  P  and  a  deviator  component  s,  where  -P  is  the  mean  of 
three  stresses:  -P  =  (1/3)*  (oj  +  4  <r^). 

*1  =  "  P  +  3 1  b-j  =  -  P  +  s  j 

-  P  +  also  —  -  P  +  s ^  (2) 

cr3  =  -  P  +  s3  cr3  =  -  P  +  S3. 


The  usual  notation  is  followed  here  such  that  stresses  are  >0  in  tension  and 
<  0  in  compression,  which  is  just  the  opposite  for  pressure;  hence  the 
negative  sign  in  front  of  the  pressure.  We  will  define  the  mean  normal  strain 
as : 

0  =  (6j  +  6 2  +  €  )  also  0  =  (€-  +  +  €3)o  (3) 


The  normal  components  of  the  strain  deviators  are  defined  as: 


01  =  ‘]  *e 
02  "  62  *  9 
‘  'i  ‘  ® 


also 


61  =  <1  -9 

92  ~  *2  "  9 
9  3  =  *3  -  0. 


(4) 


From  the  equation  of  continuity  we  have: 

.  .  .  .  V 

€1  +  e2  +  e3  =  V  • 


(5) 


It  follows  that: 

h  + W° 


and 


0-II 

9  3  V 


Using  the  above  definitions  we  may  now  write  Hooke's  law  [Eq„  (T)]  as: 


si  =  2jA (ei  “  3  v) 
s2  =  2hL^e 2  "  3  ?) 

S 3  =  Z[±  (^3  "  3v) 

p  =  _  K  — 
iV  y 


=  bulk  modulus. 


(6) 
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Also  from  equations  5  and  6,  it  follows  that: 
»1  +  “2  +  's3  =  0 

and 


s 


1 


+  s2  +  s3  =  ° 


(7) 

(8) 


which  says  that  the  distortion  components  of  the  stresses  do  not  contribute  to 
the  average  pressure. 


Bo  Plastic  Flow  Region 

The  yield  condition  of  R.  Von  Mises  is  used  to  describe  the  elastic 
limit  (see  Ref.  4  for  an  English  translation  of  the  Von  Mises  paper).  When 
the  principal  stresses  are  known,  the  yield  condition  can  be  written  as: 

(o-j  -  er2)2  +  (c -  a 3)2  +  (c -  o^)2  -  2(Y°)2  (9) 

where  is  the  yield  strength  in  simple  tension. 

The  left  side  of  this  expression  is  proportional  to  the  elastic  energy  of 
distortion  per  unit  volume  or  the  energy  required  to  change  shape  as  opposed 
to  the  energy  that  causes  a  volume  change  (Ref.  5,  p.  210).  The  expression 
states,  therefore,  that  plastic  flow  begins  when  the  elastic  distortion  energy 
reaches  a  limiting  value  [(Y^)2/ 6p]  and  that  this  energy  remains  constant 
during  the  plastic  flow.  Thus,  by  the  term  "elastic -plastic"  is  meant  the 
state  whereby  the  distortion  (change  in  shape)  component  of  the  strained 
material  has  been  loaded,  following  Hooke's  law,  up  to  a  state  where  the 
material  can  no  longer  store  elastic  energy.  All  subsequent  distortion  will 
produce  plastic  flow  and  plastic  work  will  be  done.  The  left  side  of  Eq.  (9) 
can  also  be  interpreted  in  terms  of  shear  strength.  There  are  several  ways 
of  viewing  Eq.  (9),  but  the  point  here  is  that  at  the  elastic  limit  the  left  side 
is  equal  to  a  constant.  We  have  chosen  to  interpret  the  constant  in  terms  of 
the  yield  strength  in  simple  tension  Y  If  the  tension  is  applied  in  the  <r j 

direction  and  the  lateral  stresses  au  and  c are  zero,  then  Eq.  (9)  gives 

0  ^  J 
cr ^  =  Y  .  The  simple  tension  implies  two-dimensional  flow  since  in  order 

for  the  lateral  stresses  to  be  zero  there  must  be  strains  in  the  lateral 

directions.  In  fact  the  ratio  c  / for  this  case  is  Poisson's  ratio,  also,  it 

is  noted  that  Equ  (9)  implies  that  the  yield  strength  in  tension  is  the  same  as 

the  yield  strength  in  compression  (absence  of  Bauschinger  effect). 
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In  orT  ,  or_  ,  or_  space,  Eq.  (9)  describes  the  surface  of  a  straight 
circular  cylinder.  The  axis  of  the  cylinder  is  equally  inclined  to  the  er ^ 

(r  system  of  coordinates  as  shown  in  Fig.  la.  The  radius  of  the  cylinder  is 

Y  Q 

72/3  Y  .  We  are  going  to  use  the  principal  stress  deviators  such  that 
sl+s2+s3  =  0[see  Eq.  (8)]  .  This  is  the  equation  of  a  plane  through  the 
origin  of  the  axes  of  the  principal  stresses.  The  intersection  of  this  plane 
with  the  cylinder  of  Eq.  (9)  results  in  a  circle  (see  Fig.  la).  If  the  stress 
deviators  s  j  ,  s^,  s^  give  a  point  inside  the  circle,  the  material  is  within  the 
elastic  limit. 

When  the  material  is  loaded  beyond  the  yield  strength  and  subsequently 
unloaded,  only  the  elastic  distortion  energy  is  recovered.  The  work  done 
against  the  material  while  in  the  plastic  state  is  not  recovered.  Another  way 
of  stating  this  is  that  the  loading  and  unloading  paths  are  not  the  same  when 
the  material  has  been  loaded  beyond  the  elastic  limit  (in  Fig.  2a  the  loading 
path  is  OA.B,  the  unloading  path  is  BC).  It  has  been  demonstrated  by  D.  C. 
Drucker^  that:  The  work  done  on  the  material  during  a  loading  and  unloading 
cycle  must  be  positive  or  zero;  zero  only  when  purely  elastic  changes  take 
place.  Furthermore,  the  plastic  strain  increment  must  be  normal  to  the 
yield  surface  that  separates  the  elastic  and  the  elastic -plastic  states.  We 
will  describe  plastic  flow  by  maintaining  the  stress  deviators  (s^,  s^,  s^)  at 
the  elastic  limit.  In  Fig.  lb  the  stresses  are  shown  at  state  n  and  after  an 
incremental  strain  we  will  consider  that  the  stresses  have  changed  to  state 
n  +  1.  However,  state  n  +  1  is  outside  the  yield  circle  and  our  assumption  is 
that  this  state  can  not  be  reached.  Instead,  we  will  consider  that  the 
material  flows  plastically,  but  the  stresses  remain  at  the  elastic  limit  on  the 
yield  circle.  The  plastic  component  of  strain  is  perpendicular  to  the  yield 
curve  and  it  is  the  stress  associated  with  this  component  of  strain  that  we 
want  to  limit.  Therefore,  the  new  stress  state,  instead  of  being  n  +  1  ,  is  the 
point  that  is  reached  by  a  vector  from  n  +  1  and  perpendicular  to  the  yield 
circle.  The  one-dimensional  analogy  is  seen  in  Fig.  2b  where  the  stress, 

-s  j  ,  has  a  maximum  value  for  all  strains  beyond  the  elastic  limit  point  A. 
Thus,  to  summarize  the  yield  assumption: 

(Sj  -  s2)2  +  (s2  -  s3)2  +  (s3  -  Sl)2  <  2(Y°)2 

sl  +  s2  +  s3  =  ° 


(10) 
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Fig.  2.  One-dimensional  strain  for  a  perfectly  plastic  material. 


UCRL-7322 


-8- 


which  can  be  written  as: 

S1  +  S2  +  S3  -  2/^3  (Y°)2  (11) 

If  an  incremental  change  in  the  stresses  in  an  element  results  in  a 
violation  of  the  inequality,  then  each  of  the  principal  stress  deviators 
(Sj,  s ^ 5  s3)  must  be  adjusted  such  that  Eq„  (11)  is  again  satisfied.  Hooke's 
law  is  used  to  calculate  the  stress  deviators  [Eq.  (6)].  If  a  point  falls  out¬ 
side  the  yield  circle  (Fig.  lb)  it  is  brought  back  to  the  circle  along  the  radius 
vector  of  the  point  and  hence  perpendicular  to  the  yield  circle.  This  is 
accomplished  by  multiplying  each  of  the  stress  deviators  (Sj,  s^,  s3)  by 
Jz/3  Y°  yjs 2  +  s ^  +  s2.  By  adjusting  the  stresses  perpendicular  to  the  yield 
circle,  only  the  plastic  components  of  the  stresses  are  affected.  The  observed 
incompressibility  of  the  plastic  state  is  implicit  in  this  procedure.  Note  that 
there  is  always  a  background  pressure  stress  present,  whether  the  material 
is  in  an  elastic  or  an  elastic -plastic  state,  but  it  is  independent  of  the  plastic 
flow.  This  is  in  agreement  with  the  observed  behavior  of  ductile  metals. 

The  above  formulation  corresponds  to  a  perfectly  plastic  material,  i.  e.  , 
material  that  flows  plastically  under  a  constant  stress  without  work-hardening 
(see  Fig.  2).  For  a  work-hardening  material  the  stress  (-Sj)  will  increase 
monotonically  with  strain  (Cj)  for  strains  beyond  point  A  instead  of  remaining 
constant  as  for  the  perfectly  plastic  material  shown.  Work  hardening  can  be 
introduced  into  the  calculation  by  making  the  constant  Y°  in  Eq.  (9)  a  function 
of  the  strain  energy,  for  example.  Also,  when  enough  work  has  been  done  to 
melt  the  material,  the  value  of  Y°  can  be  set  to  zero.  In  this  way  an  all¬ 
hydrodynamic  description  will  follow  since  the  stress  deviators  will 
automatically  be  set  to  zero  by  the  above  procedure  and  the  only  remaining 
stress  will  be  the  pressure  P.  Time-dependent  yielding  can  be  macroscopically 
represented  by  selecting  a  high  yield  constant  Y  if  the  strain  rates  (€j  >  e2’ 
e„)  are  above  some  prescribed  value. 

3  .  o 

In  the  negative  pressure  region,  the  pressure  is  cut  off  at  P  =  -  (1/3)Y 

consistent  with  a.  simple  tension  test. 
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The  complete  equation  of  state  is  given  by: 


°‘l  =  -P+Sl 


(i) 


si  =  2jx[ei  "  Tv] 
ff2  =  “  p  +  s2  (ii)  {  s2  =  2^[€2  "  3  v] 

s3  =  2p[€3  - 


r3=  “P  +  S3 


(iii) 

2,2,2 
S1  +  S2  +  S3 

<  |  (Y°)2 

(iv) 

S1  +  s2+  s  3 

=  0 

(v) 

P  =  P(V) 

(12) 

_  1  Yo 

"3 

(vi) 

Minimum  P 

C.  Experimental  Equation  of  State 

Consider  a  one-dimensional  shock  traversing  a  material  such  that  there 
is  a  strain  in  the  X  direction  and  zero  strains  in  the  Y  and  Z  directions. 

This  is  the  geometry  whereby  Hugoniot  equation- of- state  data  are  obtained. 

A  shock  exists  that  takes  the  material  through  an  elastic  to  an  elastic -plastic 
stated 

For  one -dimensional  flow,  the  X,  Y,  Z  coordinates  are  the  principal 


The  total  stress  in  the  X  direction  is: 

0-  =  -  P  +  s  .  U4) 

X  X 

We  will  assume  it  is  cr^  that  is  obtained  by  Hugoniot  measurements. 


X-X  -  Hugoniot 

A  -  Hugoniot  elastic  limit 

9  =  p/p°  =  VV 

pO  =  reference  density 
p  =  actual  density 


Fig.  3 
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For  1-D  flow  the  equation  of  continuity  gives  =  V/V.  The  complete 
equation  of  state  for  one-dimensional  geometry  is  described  by: 


(iv)  s  +  2s  =0 

x  y 

(v)  P  -  P(V) 

(vi)  Minimum  P  “  -  ~  Y^.  (15) 


Up  to  the  elastic  limit,  point  A,  we  have: 


v 

p  =  -  K— 

jr 

and  s  =  2|i  [V/V  -(1/3)  V/V] 

X 

P  =  -  K  In  V 

and  s  =(4/3)p  In  V  =  -  2  s 

x  y 

0-  =  KlnV  +  |plnV 

x  3 

/  \ 

=  K  +  V. 

At  point  A  : 


(16) 


or 


(24/9) (P  In  VA)2  =  |  (Y°)2 


Y°  =  2 p  In  V, 


2(i(o-x)A 


(17) 


X  +  2p 

This  gives  the  maximum  yield  strength  Y  if  the  Lame  constants  and  the 
Hugoniot  elastic  limit  are  known. 

For  points  beyond  A  we  have:  (s^  +  2s^  =  (Y  ) 
which  reduces  to: 


s  =  ±  \  [from  (iv)  Eq.  (15)]  9  (18) 

x  3 

[if  the  material  yields  in  tension,  s  >0,  in 

0  x 

compression,  sx  <  0  (Y  is  always  >0)] 


-11- 


UCRL-7322 


so  the  total  stress  resulting  from  a  shock  from  =  0  to  a  value  above  point 
A  is: 

-  «r  =  P(V)+-|(Y°)  (19) 

x  3 

P(V)  can  be  expressed  conveniently  as: 

P(ti)  =  a(t]  -  1)  +  b(q  -  1)^  +  c(t)  -  1)^.  (20) 

Here  a,  b,  and  c  are  constants  such  that  P(q)  +  (2/3)  Y^  reproduces  the 
Hugoniot  for  shocks  above  the  elastic  limit  and  P(tj)  =  -Kin  V  for  pressures 
below  the  elastic  limit. 


Fig.  4 


The  result  of  using  an  equation  of  state  as  given  by  Eqs.  (15)  is  a  loading 
path  O,  A,  B  (Fig.  4)  and  an  unloading  path  B,  C,  D.  Experiments  on  metals 
in  the  low-pressure  range  (0  -►  50  kb)  have  demonstrated  the  difference 
between  the  hydrostatic  P(r|)  and  the  Hugoniot  (cr  )  curves.  At  high  pressures 
(hundreds  of  kilobars)  for  some  metals  the  sound  speed  behind  the  shock  has 
been  measured  to  be  of  the  order  of  20%  faster  than  that  predicted  by  hydro- 
dynamic  theory.  This  gives  reason  to  extend  the  low-pressure  model  up  to 
high  pressures.  Upon  unloading  from  a  high  pressure,  the  material  unloads 
first  elastically  along  BC;  the  slope  of  this  path  is  characteristic  of  the 
elastic  unloading  velocity.  Consequently,  the  rarefaction  travels  faster  than 
would  be  the  case  if  the  material  unloaded  entirely  along  the  P(t])  path. 
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Experiments  that  have  measured  high  sound  speed  behind  shocks  in  metals 

g 

have  been  performed  in  Russia,  Stanford  Research  Institute,  and  at  this 
Laboratory. 

9 

Using  Hugoniot  data  of  J.  M.  Walsh  and  the  elastic  data  of  C.  D. 
Lundergan,^  the  constants  for  Eqs.  ( 1  2)  for  aluminum  are: 

Y°  =  0.002976  mb  (from  Eq.  (17)  with  (<rJA 

p  =  0.248  mb  =  0.0063  mb  and  =  0.994) 

P  =  0.73  (n  -  1)  +  1*72  (t)  -  l)2  +  0.40  (^  -  l)3 

0  ->7 

p  =  2o  7. 

Figures  5  and  6  show  the  results  of  finite  difference  calculations  of  a 
flying  aluminum  plate  striking  a  target  plate-  The  calculation  shows  that 
even  though  the  yield  strength  is  small  compared  to  the  total  stress  c r  , 
the  effect  on  the  wave  is  very  pronounced- 

In  principle,  calculations  of  this  type  in  conjunction  with  experiments 
could  be  used  in  determining  the  properties  of  materials  at  high  pressure 
based  on  the  model  described  by  Eqs.  (12).  Front- surface  velocity  measure¬ 
ments  for  various  thicknesses  of  target  plates  could  determine  when  the 
elastic  wave  overtakes  the  shock  front  (see  Fig.  6)  and  this  would  establish 
the  slope  of  BC  in  Fig.  4.  The  magnitude  of  the  step  behind  the  shock  in 
Fig.  6  would  correspond  to  point  C  of  Fig.  4. 

If  the  material  behaves  entirely  hydrodynamically ,  then  the  step  behind 
the  shock  front  will  not  be  present.  In  Fig.  5  it  is  seen  that  the  rarefaction 
due  to  the  hydrodynamic  unloading  is  proceeding  much  slower  than  the  elastic 
unloading.  The  hydrodynamic  rarefaction  however,  can  be  seen  to  be  over¬ 
taking  the  shock  front  from  the  increase  in  the  slant  of  the  rear  of  the  wave 
as  time  increases. 

It  should  be  noted  that  to  account  for  high  sound  speeds  behind  a  shock 

it  is  only  necessary  to  postulate  that  the  material  unloads  first  elastically 

and  then  plastically.  The  result  would  still  follow  even  though  at  high 

pressures  the  -  <r  curve  merged  with  the  P(t])  curve.  We  have  maintained 
0  ^ 

the  values  of  Y  and  p  constant  in  these  calculations  since  the  details  of  the 
elastic  unloading  from  high  pressures  are  not  known. 
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5.  Stresses  from  a  flying  plate  striking  a  target  plate. 
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PART  II.  ONE -DIMENSIONAL  ELASTIC -PLASTIC  FLOW 

For  time -dependent  flow  in  one  space  variable  (r) ,  the  principal 
equations  for  plane  (d  =  1),  cylindrical  (d  =  2),  and  spherical  (d  =  3) 
geometries  are: 


equation  of  motion 

o*  as 

P  £  = _ £ 

V  9r 


+  (d  -  1) 


Sr  =  -(P+  q)  +  Sj 
E0=  -(P  +  q)  +  s2 


(21) 


equation  of  continuity 

V  _  1  9(rd~ 1 U) 

V  d-1  9  r 

r 


(22) 


energy  equation 

E  -  V 


S1  M  +  (d  -  ])  s2  k2 


artificial  viscosity  (linear  q) 
^0  p°  /9U\  . 

* =  c  V  Ira  Ar 


+  (P  +  q)  V  =  0 


=  constant 


(23) 


(24) 


equation  of  state 


anisotropic 

stresses 


(25) 


NOTE:  Three  stresses  are  identified  here,  even  though  they  are 
not  all  required,  so  as  to  maintain  an  analogy  with  the 
two-dimensional  calculations  in  Part  III. 


velocity 

strains 


f  .  _  9U 

€1  3r 

U 

=  — 

2  r 

^  €3  =  e  ^  for  d  =  3 

€3  =  0  for  d  -  2 

=  0  for  d  =  1 

V  3  z 
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hydrostatic  pressure 


P  =  a(ri  -  1)  +  b(t|  -  l)2  +  c(r)  -  l)3  +  d-^E 


■n  = 


v 


_p_ 

o 

p 


Von  Mises  yield  condition 


(s2  +  S2  +  s2)n+1  -  (2/3)(Y°)2  <0  Y°  =  material  strength 


Notation: 


r 

U 

zr’  2e 
sr  V  s 3 


€r  V  e3 


v 


p 

p 

E 


0 


space  coordinate 
velocity  in  r  direction 
total  stresses 
stress  deviators 


strains 

hydrostatic  pressure 
relative  volume 
reference  density 
actual  density 

internal  energy  per  original  volume 
The  dot  over  a  parameter  signifies  a  time  derivative  along  the 
particle  path. 


Application 

The  finite  difference  equations  for  the  above  set  of  partial  differential 
equations  are  given  in  Appendix  A. 

Figure  5  shows  a  calculation  in  plane  geometry  (d  =  1)  of  a  flying  plate 
coming  from  the  left  with  a  velocity  U  =  0*08  cm/ps ec  and  striking  a  target  at 
rest.  The  materials  are  aluminum,  with  the  constants  given  in  Part  I.  At 
this  impact  velocity  the  elastic  signal  travels  faster  than  the  shock  speed  and 
the  shock  breaks  into  two  components.  In  Fig.  6  the  same  problem  has  been 
calculated  with  a  higher  impact  velocity  (0.2  cm/psec).  In  this  case  the  shock 
speed  is  greater  than  the  elastic  signal  speed  and  the  shock  front  is  the  same 
as  though  the  material  were  following  an  all  hydrodynamic  descriptions. 
However,  it  can  be  seen  that  the  elastic  relief  wave  that  originated  from  the 
rear  of  the  flying  plate  is  overtaking  the  shock  front. 
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Figure  7  shows  the  stress  waves  resulting  from  a  1-cm-radius 

spherical  charge  of  high  explosive  (comp.  B)  detonated  at  the  center.  The 

region  from  1  cm  to  5  cm  is  treated  as  an  elastic  material  with  P  =  K(ri  -  1), 

K  =  1.39  mb,  p°  =  8.9,  pi  =  0.46  mb,  and  Y°  =  oa.  The  stress  waves  are 

traveling  at  a  velocity  C  =  [(K  +  -|  p)/ p*  V/2  cm/p  sec.  A  second  shock  can  be 

J  .12 

seen  that  has  originated  from  the  high- explosive  cavity. 

The  second  shock  is  a  hydrodynamic  effect  and  is  not  a  result  of  the 
elastic  property  of  the  material.  The  radial  stress  is  seen  to  be  followed  by 
a  tension  tail  of  about  -10  kb.  The  tension  is  due  entirely  to  the  elastic 
property  of  the  material.  If  the  material  were  described  by  hydrodynamics 
alone ,  there  would  not  be  a  tension  portion  behind  the  outward-travelling 
pressure  stress  wave. 


LOBARS 
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PART  III.  TWO-DIMENSIONAL  ELASTIC -PLASTIC  FLOW 


The  equations  of  motion  listed  below  are  those  used  by  the  HEMP  code, 

a  program  that  solves  the  equations  by  finite  difference  techniques  on  the 

IBM  7030  electronic  computer.  The  derivation  of  the  equations  can  be  found 

in  Ref.  13.  The  problem  is  formulated  in  Lagrange  coordinates  with  sliding 

interfaces  allowed  between  an  elastic  and  a  hydrodynamic  region  or  between 

two  hydrodynamic  regions,  but  not  between  two  elastic  regions.  However, 

an  elastic  region  may  slide  along  a  fixed  boundary.  The  equation  of  state  is 

used  in  the  same  manner  as  described  in  the  preceeding  sections.  There  is, 

however,  the  additional  complication  that  the  stress -strain  relationship  must 

be  independent  of  a  rigid  motion  and  hence  the  incremental  stress-strain 

relationship  must  be  corrected  for  a  rotation  in  the  x-y  coordinate  system 

(Ref.  14).  When  a  zone  is  displaced  from  an  initial  state  of  stress,  there 

may  be  a  rotation  through  an  angle  co  as  well  as  a  distortion.  The  rotation 

will  not  contribute  to  an  increase  in  stress,  but  the  state  of  stress 

(s  ,  sn  ,  Tn  )  originally  in  the  zone  has  been  rotated  through  the  angle  w. 

\xx  yy  xy / 

Since  the  equations  of  motion  are  referred  to  the  fixed  x-y  coordinate  system, 

the  rotated  stresses  must  be  recalculated  in  terms  of  the  coordinate  system. 

The  transformation  equations  (Ref.  5,  p.  110)  result  in  a  correction  6  that  is 

added  to  the  stresses.  The  stresses  can  then  be  incremented  by  the  strain 

, n  ...  , n+ 1  .  .  .  ..  . n+1 


that  occurred  between  time  t  and  time  t 
The  rotation  angle  is  given  by: 

=  I  _  ^l\  Atn+1//z 


to  give  the  stresses  at  time  t 


sm  co 


Z  \3x  3y i 


It  is  not  practical  to  increment  the  stresses  in  the  principal  stress 
coordinate  system  because  the  principal  stress  directions  are  not  unique  and 
also  there  can  be  large  changes  in  directions  between  two  consecutive  time 
steps.  Therefore,  the  operation  of  transforming  to  the  principal  stress 
coordinate  system  and  then  back  to  the  x-y  coordinate  system  every  cycle 
would  become  complicated  and  inaccurate.  Even  though  the  program  to  be 
described  does  make  the  transformation  to  the  principal  stress  system  to 
test  for  yielding,  the  directions  of  the  principal  stresses  are  not  required  in 
the  calculations  and  the  complications  involving  these  directions  is  avoided. 
(The  Von  Mises  yield  condition  can  be  used  without  reference  to  the  principal 
stresses  (see  Appendix  B  7-b),  however,  the  geometrical  interpretation  of  the 
yielding  is  easier  if  the  principal  stress  deviators  are  used  instead  of  the 
stress  deviators.  ) 
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A.  Basic  Equations  in  the  HEMP  Code 

Equations  of  motion  in  x-y  coordinates  with  cylindrical  symmetry  about 
the  x-axis  (it  is  desirable  to  have  the  problem  formulated  also  in  plane  x-y 
coordinates;  for  this  case  the  terms  marked  |  ‘  are  set  =  0): 


92  9T 

xx  ,  _ x; 

9x  9y 


T 

+  -S 

y 


=  px 


z  -  Z, 


Z  =  S  -  (P  +  q), 

XX  XX 


Z  =  s  -  (P  +  q) » 

yy  yy 

=  see  ‘  (p  +  <i)' 


Equation  of  continuity: 

1  =  +  x 

V  9x  9y  y 


Energy  equation: 


E=-(P+q)V  +  V(s  €  +  s  k  +s  e  +T  k 

'  H/  V  xx  xx  yy  yy  oo  do  xy 


xy  xy 


Artificial  viscosity:  (quadratic  "q") 
q  =  C  2  p°(Y/V)2  A/V 


where 


C  =  constant 
A  -  zone  area 

=  reference  density. 


Equation  of  state: 


stress 

components 


0 

yxx  “  3  V  J 


&  =  2p' 


"syy  =  2p(Vy  -  iv) 

see  =  ~  sv) 


T  =  p(e  )  +  6  , 

xy  xy  xy 
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where 

|x  =  shear  modulus 

6  =  correction  for  rotation  (see  text). 


velocity 

f  .  _  9x  . 

I  Sex  9x  €00 

_  y 
y 

strains 

|  j  =^L  k 

=  9x  +  9x 

L  yy  9y  xy 

9x  9y 

Hydrostatic  pressure 

P  =  a(ri  -  1)  +  b{r\  -  l)2  +  c(t]  -  l)3  +  dr|E, 
T|  =  l/V  =  p/p°. 


Von  Mises  yield  condition 

(s  2  +  s  2  +  S2)  -  (2/3)(Y0)2  <  o 

where 

=  material  strength 

(s  ,  s2,  s3)  are  the  principal  stress  deviators. 


Notation: 


x,y  space  coordinates 

x  velocity  in  x  direction 

y  velocity  in  y  direction 

22  ,  21  ^  total  stresses 

xx  yy  uu 

T  shear  stress 

xy 


s  ,  s 
xx  yy 

€  ,6 

xx  yy 

P 

V 

E 

P 


3  sQn  stress  deviators 

,  €  strains 
00  xy 

hydrostatic  pressure 
relative  volume 

internal  energy  per  original  volume 
density 


(31) 


(32) 


(33) 


The  dot  over  a  parameter  signifies  a  time  derivative  along  the  particle  path. 
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B.  Finite  Difference  Scheme 


The  following  integral  definitions  of  the  partial  derivatives  are  used 
(seep.  327,  Ref.  15): 

f  F(n  •  i)  dS 

9 E  =  Z£ _  (34) 

3x  lim  A 
A~^0 


C  F$  •  JA)  dS 

3F  _  _^C _ _ 

3y  lim  A. 

A—0 

where 

C  is  the  boundary  of  area  A 
S  =  arc  length 
n  =  normal  vector 
t  =  tangent  vector. 


(35) 


A  _  9x  A 

n  3n  1 


(36) 


Applying  the  above  to  the  quadrilateral  1, 
F  defined  at  the  points  I,  2,  3,  and  4: 


2,  3,  4,  we  have  for  a  parameter 


A  =  area  of  quadrilaterial 


Cf({.  •  i)dS  =  +  Jr|gdS  (37) 

=  -  [F23(y 2  '  y3>  +  F34(y3  '  y4>  +  F41(y4  '  yl>  +  F12(yl  ‘  y2>]  (38) 

where  =  (F2  +  F3)/2  5  etc. 

S°  -  X  [F23(y2  -  y3>  +  F34(y3  '  y4>  +  F41(y4  '  yI>  +  F12(yl  ‘  y2>]  <39) 

=  +  -k  K  -  F4»y3  “  yl>  -  <y2  -  y4)(F3  -  Fl>]  • 
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Similarly 

If  =  '  TA  (F2  "  F4)(x3  "  Xl)  "  (x2  '  X4)(F3  '  Fl}  •  (40) 

The  quantities  3F/9x  and  9F/3y  are  considered  to  be  defined  at  the  center  of 
the  quadrilateral. 

Using  the  above  difference  equations  we  can  now  write  an  expression 
for  3x/3x  and  9y/9y  at  a  given  time  and  position.  The  difference  scheme 
to  be  used  defines  the  velocities  (x  and  y)  at  l/z  time  increments  and  the 
space  quantities  (x  and  y)  at  integral  time  increments. 

If  we  use  the  definitions: 

xn+l/2  =  ,/2  (xn+l  +  xnK  etc 

and 

A”+l/2  =  1/2  (An+1  +  An) 

where  An  and  A.n+^  =  area  of  the  quadrilateral  at  time  tn  and  time  tn+*  , 
respectively,  then  the  difference  equations  will  give: 

(continuity  equation  in  plane  x-y  (41) 

geometry  where  A./ A  =  V/V). 

It  is  obviously  very  desirable  that  a  difference  scheme  have  this 
property  since  it  leads  to  zero  truncation  error  in  the  numerical  integration 
of  each  of  the  terms  in  Eq0  (41). 

We  will  now  consider  the  continuity  equation  in  x,  y  coordinates  with 
cylindrical  symmetry  about  the  x  axis. 


3x  9y  ,  A. 

97  exactly  A 


+  i  =  i . 

3x  9y  y  V 


(42) 


Here  V  is  the  volume  swept  out  when  the  area  A  is  rotated  about  the  x  axis. 


V  =  yaAa  +  VS. 

A.  =  area  of  Aa  A.  =  A.  +  A., 
a  a  b 


A,  =  area  of  Ab 
b 

ya  =  1/3  (y2  +  y3  +  y4} 

yb  =  1/3  (y2  +  y4  +  y2) 
ya  =  i/3  (y2  +  y3  + 

yb  =  1/3  (y2  +  y4  +  yj) 


(43) 
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A  very  good  approximation  for  the  third  term,  y/y,  of  Eq.  (42)  is  given 


by: 


X  = 


A  n+ 1 / 2 •  ,  .n+l/2. 

Aa  ya  +  Ab  ^b 


y  An+l/2_n+l/2  ~  A“+1/2yf+1/2 


(44) 


It  is  important  to  recognize  that  the  difference  equations  for  the  terms 
in  brackets  in  Eq.  (42)  are  the  same  as  for  the  left  side  of  Eq.  (41)  i.  e.  , 
they  are  independent  of  the  coordinate  system.  Since  there  is  essentially 
zero  truncation  error  with  the  integration  of  these  terms,  it  is  possible  to 
calculate  V/V  from  the  coordinates  and  express  y/y  as: 


1  =  1 
y  V 


dx  dy 
dx  dy 


(45) 


For  the  acceleration  routines,  the  given  parameter  F  in  Eqs.  (34)  and 
(35)  is  defined  at  the  center  of  a  quadrilateral.  The  area  enclosed  in  the 
integration  is  the  area  I,  II,  III,  IV  in  the  figure  below. 


The  corresponding  difference  equations  for  the  l  and  f  components  of 
acceleration  become: 

Cf(&  .i)dS=.  [j?Q(yu  -  ym>  +  F@(ym  -  yIV)  +  F@(yIV  -  Yj)  +  -  yn)]  . 

(46) 

^F(n  •  j)  dS  =  +  [F0(xI];  -  Xjjj)  +  F0(xn][  -  x^)  +  Fq{x1v  -  Xj)  +  F^Xj  -  xn)J  , 

(47) 

The  area  I,  II,  III,  IV  is  considered  to  be  the  mean  of  the  quadrilateral 
areas  Aq,  A^,  .A^,  and  A^.  As  can  be  seen  in  Appendix  B,  the  quadrilaterals 
are  weighted  by  the  four  corresponding  densities. 
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C.  Applications 

The  complete  set  of  difference  equations,  including  sliding  interfaces 
and  boundary  conditions,  is  given  in  Appendix  B.  In  the  following  figures, 
applications  of  these  equations  to  specific  problems  are  shown.  The  plots 
were  made  directly  from  the  high-speed  computer  by  a  cathode  ray  tube  and 
then  photographed.  Zones  are  shaded  with  an  intensity  weighted  by  the  rate 
of  zone  compression.  Hence,  the  shocks  or  detonation  fronts  are  traced  on 
the  grid. 

Figure  8  shows  the  time  sequence  of  events  due  to  a  charge  of  explosive 
detonated  at  constant  volume  against  a  plate  of  copper.  The  horizontal  line 
through  the  middle  of  the  grid  is  an  axis  of  cylindrical  symmetry.  The 
equation  of  state  for  the  copper  was  derived  in  the  same  manner  as  previously 
discussed  for  aluminum.  The  yield  strength  used  was  Y  =10  kb  and  shear 
modulus  |jl  =  460  kb. 

Figure  9  shows  the  directions  of  the  maximum  principal  stresses  for 
the  above  problem.  The  maximum  principal  stress  with  the  sign  convention 
being  used  is  the  tension  component  of  the  stresses.  Therefore,  the  line 
segments  point  in  a  direction  normal  to  the  direction  of  an  incident  shock. 

By  examining  the  plots,  the  progress  of  rarefactions  from  the  free  surfaces 
can  be  followed  since  the  lines  point  in  the  direction  of  tensions. 

Figure  10  shows  the  same  problem,  but  the  copper  is  described  by  a 
hydrodynamic  equation  of  state  with  no  strength  of  materials.  It  can  be  seen 
that  the  crater  is  much  deeper  than  in  the  preceding  problem.  The  crater  lip 
is  seen  to  increase  with  time.  With  a  hydrodynamic  description  there  are  no 
restoring  forces  if  the  material  deforms  at  constant  volume.  The  elastic- 
plastic  description  gives  rise  to  restoring  forces  that  resist  a  change  in  shape 
even  though  the  volume  remains  constant,  and  hence  no  lip  is  formed  for  the 
problem  shown  in  Fig.  8.  A  lip  will  form,  however,  for  calculations  made 
with  yield  strengths  of  1  to  2  kilobars. 

Figures  11  and  12  show  the  time  sequence  of  stresses  in  a  copper  plate 
resulting  from  the  interactions  of  detonation  fronts  in  a.  high  explosive  (PBX 
9404).  The  calculation  was  made  in  plane  geometry  and  the  detonation  centers 
are  lines  perpendicular  to  the  page  on  the  right-  and  left-hand  lower  corners. 
The  left-  and  right-hand  boundaries  are  planes  of  symmetry.  The  copper  has 
a  yield  strength  =  10  kb  in  this  calculation. 
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t  =  0.5//sec 
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t=  l.5//sec  t  =  2.25//sec 

GLL  -639-2185 

Fig.  11.  Stresses  induced  in  a  copper  plate  by  high-explosive  detonations. t 
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TRANSMITTED  SHOCK 


t=  2.99//sec 


t=  4.03//sec 


t  =  4.47//sec 


GLL -639-2186 


Fig.  12.  Continued  from  Fig.  11. 
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The  interaction  of  the  two  detonation  fronts  results  in  the  trailing 
shocks  with  high  localized  pressures  in  the  region  of  the  cusp.  The 
detonation  front  can  be  seen  to  transmit  a  shock  into  the  copper  and  reflect  a 
shock  into  the  high  explosive.  The  trailing  shocks  induce  a  circular  shock 
that  sweeps  from  the  left  and  right  of  the  line  of  interaction.  This  shock  can 
be  seen  to  punch  out  the  center  of  the  copper.  In  the  last  frames  the  reflection 
of  the  trailing  shock  from  the  fixed  boundaries  can  be  seen.  The  dark  shading 
where  the  reflected  shock  from  the  copper  meets  the  trailing  shocks  indicate 
a  high-pressure  region. 

Figure  13  shows  the  directions  of  the  maximum  principal  stresses  in 
the  copper  for  the  above  problem.  The  expanding  circular  shock  due  to  the 
trailing  shocks  can  be  seen  to  be  progressing  through  the  already  shocked 
copper. 

Figure  14  shows  a  comparison  of  the  above  problem  with  a  calculation 
made  with  no  material  strength  in  the  copper.  This  calculation  indicates  a 
low  density  region  along  the  front  surface  as  well  as  in  the  large  circular 
section  in  the  middle. 

Figure  15  shows  the  explosion  of  an  iron  cylinder  that  had  a  charge  of 
PBX  9404  inside  that  was  detonated  from  the  right-hand  surface.  On  the  left, 
the  iron  cylinder  motion  has  been  calculated  with  no  strength  of  material  and 
on  the  right,  the  motion  was  calculated  with  a  yield  strength  =10  kb. 
Comparing  the  two  calculations,  it  can  be  seen  that  the  elastic -plastic  version 
shows  the  thickness  of  the  iron  to  be  slightly  dilated.  Also,  the  end  of  the 
cylinder  has  maintained  more  of  its  original  shape  compared  to  the  hydro- 
dynamic  version. 
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This  work  was  performed  under  the  auspices  of  the  U.  S.  Atomic 
Energy  Commission. 
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FIXED  BOUNDARY  > 
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HE.  DETONATION 
SURFACE 


AXIS  OF  CYLINDRICAL  SYMMETRY 


HCMP —  CYCLE  TIME 

+52  90000  +5 1  •  20000  +4-3  ■  99999 


DETONATION  FRONT 


HEMP  — 

+52  90000 


CYCLE 
+53 ■ 22 1 00 


TINE 

+52  I  1 04-9 


HEMP  — 

+52  90050 


CYCLE 
+53 • 22 1 00 


TINE 

+52  I  1 04:9 


HEMP - -  CYCLE  TINE  HEMP--  CYCLE  TINE 

+52  90000  +53  23300  +52  1 4-393  +52  -  90050  +53  ■  29300  +52  •  \  4-399 

GLL-639-2189 

Fig.  15.  High -explosive  detonated  inside  of  an  iron  cylinder,  {left:  Iron 
treated  as  a  hydrodynamic  material;  right:  Iron  treated  as  an  elastic-plastic 
material  with  a  yield  strength  Y®  =  0.010  mb). 
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APPENDIX  A 

FINITE  DIFFERENCE  EQUATIONS  FOR  THE  EQUATIONS 

GIVEN  IN  PART  II 


The  material  is  divided  into  mass  intervals 

,*(&- Ml 


mj+l/2 


plane:  d  =  1 

cylinderical:  d  ~  2 
spherical:  d  =  3 


j  =  1,2,.  .  .N  (see  Ref .  1 6). 


1.  Equation  of  motion 

(a)  U?+1/2  =  U?"1/2  + 

'  7  -l 


i  n 


i  n 


lj+l/2 


r  n 


j-  1/2 


+  Atn(p“j  (d  -  1) 


where: 


,  n 


■  ,—n  ,  n-  l/2v  ,  n  l 

j+l/2s  -'P+Cl  )  +  Sl  j+l/2 


,n\ 

'QJ  j+l/2 


-  (Pn  +  q 


X  n  _  1 


0 

pj+l/2 


n-l/2,  +  „n\ 
n 

rj+l  r- 


2  j+l/2 


n , 


V?H/2 


+  Pj-i/2 


,  n  n 
'r.  -  r.  , 

J.  ■  J---1- 

V"-  1/2 


KW  -  H+l/2 


1  /  n  ,  n 
r.,  ,  +  r. 
J+l  J 


L 


'.y2\ 


>°  /j+l/2 


(Zr)”-l/2  -  (Zetl/2 

1  /  n  n  N 
■s-  (r.  +  r.  , 

2  \j  j-1 


rn' 


,P0/j-l/2 


At  an  outside  regional  boundary  J 


,  n  _  0 

h  "  PJ-l/2 


n  n 

r  J  r  J  -  I  \  1 


VJ-l/2 

(Sr)j-l/2  ~  (Ze)j-l/2 

1  /'n  ,  n  \ 
2(rJ+rJ-l) 


rn\ 


f— 

\P0  Jj  - 1/2 
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At  an  inside  regional  boundary  J 


,  n  _  1  0 

4*  T  ->  p  T 


J  2  P  J+l/2  l  vn 


rJ+l  ‘  rJ 


outside 

boundary 


J+l/2 


(Zr)j+l/2  "  Fe)  J+I/2I  /Vf 


1  /  n  ,  n  N 

2  (rJ  +  rJ+l 


1 _ \  b  inside 

p0  AJ+l/2  boundary 


increasing 

j 


For  a  free  surface  at  j  =  J,  the  stresses  are  set  to  zero  at  J  +  l/2  for  an 
outside  free  surface  or  J  -  l/Z  for  an  inside  free  surface.  The  zones 
adjacent  to  an  elastic  region  are  considered  to  be  hydrodynamic. 
cu\  n+1  _  n  ,  TTn+l/2  Afn+l/2 


=  r  +  U 

J  3 


2.  Equation  of  continuity 


<a)  <1/2  = v? 


j+1/2 


n+1/2 


'j+1/2 


TTn+ 1/2  /  n+  l/z\  d-1 

j+1  j+1 


/  n  \  nt  1  __  I 

(b)  Vl/2  „n+l 


j+l/2 


-v?Hri/zY 


t  n+l/2  1  /  n+1  ,  n 

here  rj+i  =  2  h+i  +  rj+i 


x  j+1/2 


^,2k//2)3 


3.  Anisotropic  stresses 

Tn+l/2  TTn+l/2 

/  .  \n+  l/2  j+1  _ 1 _ 

€i  j+1/2  n+1/2  n+1/2 

.  '  7  rj+i  "  rj 

velocity 

StrainS  \  TTn+l/2  .  TTn+l/2 

/.  \n+ 1/2  j+  1 _ j  _ 

i  €2 )j+ 1/2  rn+I72  +  r n+1/2 
7  j+1  j 


-  0  for  d-1 


Correction  term  for  d  =  3  only 
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1 }  3+1/2 


■.);,A*4es‘ntlA 


1  /Vn+1  -  v"' 


3V  vn+1/2  /j+l/2 


stress  1/  \  .  , 

deviators 


f  V  ,  _  17-  \n+l/2  A  n+l/2 
,s2  j+l/2  +  ^  h  j+l/2  At 


1  /vn+!  -  v0' 


3V  Vn+1/2  /j+1/2 


3  j+l/2 


\n+l  ,  /  \n+l 
Sl)j+l/2  +  S 2  j+l/2 


4.  Artificial  viscosity 

n-fl/2  ^0  n+l/2  TTn+l/2  TTn+l/2i  ,  .  ,  Or0  TTn«a/c,  TTi 

qj+]/2  =  C  aP0Vl/2  j+l  "  Uj  I  calculate  only  if:  U.+  /  <  U. 


n+l/2  <TJn+l/2 

j+1  j 


C  =  constant  l/ 2 
a  =  local  sound  speed 


and:  V^/2-V"l/2)<0 


5.  Energy  equation 


dE  =  dE„  +  dE,  +  dE 

ri  1  l. 


dE  =  -  (P  +  q)  dV 


dEj  =  V  s  jl  j  At 


dE2  =  V  s2e2At(d-l) 

<*>  (=&  ■  (ei)”+./2 + <#(■,)  pJ/i  (•.)?;#  -n+1/2 


(b)  (Ez) j+l/2  ‘  (E2 ) j+ 1/2 


=  E0  “l/o  +  (d  -  1)V 


n+ 1/2  /  \  n+ 1/2  / .  \  n+ 1  /2  .  .  n+ 1/2 


1/2  S2)j+l/2  \e2 /j+l/2 


(C)  EFl/2-  - 


+  qn+1/2]  •  [Vn+ 1  _  vn]  +  dEj  +  dE2 

i  +  —  ia_-d  •  _Vn+  1  -  Vn 


j+l/2 
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This  is  a  composite  of  the  stability  criteria  given  in  the  Von  Neumann  and 
Richtmyer  paper  that  introduced  the  "qn  method  for  calculating  shocks 
(J.  Appl.  Phys.  2J_,  March  1950). 
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APPENDIX  B 


FINITE  DIFFERENCE  EQUATIONS  FOR  THE  EQUATIONS 

GIVEN  IN  PART  El 

1„  Mass  zoning  for  cylindrical  symmetry  about  the  x-axis. 

The  material  is  divided  into  quadrilaterals  with  a  grid  j-k  that  moves 
with  the  material.  In  the  figure,  the  centers  and  the  vertices  of  the  quadri¬ 


laterals  will  be  denoted  as  follows: 

1  =  j ,  k 


<D  =  j  +  L  k  +  j 


y 


®  =  j  -  J>  k  +  j 


j  "  V  k 


2 


(£)=j  +  |-,  k-j 


2  =  j  +  l,k 

3  =  j  +  l,k  +  1 


4  =  j  ,  k  +  1 


4 

. 3  | 

© 

\_a 
,b  \ 

i 

© 

2 

j- 

i  j 

jH 

hi 

k+1 


k-1 


x 


The  mass  at  time  zero  associated  with  each  quadrilateral  is  obtained  by 
the  product  of  the  initial  density  and  the  volume  swept  out  by  the  quadrilateral 
rotated  about  the  x-axes.  For  example,  the  mass  at  time  zero  for  quadri¬ 
lateral  (D  is  calculated  as  follows: 

O' 


<a)  M® = KH1 


0 


0 


y2  +y3  +y4 


0 


+  lYj0  +  y 


o 


o\  .  0 

y4  JAb 


CD 


masses  and  are  calculated  similarly. 


(b)  / 


A^  =  area  of  Aa;  A^  =  area  of  Ab: 

/  .  \  n  1  |  n  f  n  n\  ,  n/n  n\  n  j  n  n 

\  a/(p  “  2  (_x2  Y  3  -  V4  )  +  x3  VV4  ‘  y2  )  +  x4  ^2  '  ^3 

K)®  =  £[V  (v/  -  ^  (*,"  '  ^  ^  x"  (v2n  -  ^ 

A  <D  =  (Aa)fD  +  KL  ' 


Conservation  of  mass 

6 


n  1  /  p 
V®  "  3  I' M 


<D 


(n  ,  n  ,  n  \  »  n  f 

y2  +  y3  +y4  )\ 


n 


n 


+  Yi  +  +y  / 


n 


n 


CD 


v 


n 


*  l-°*n) 


©■ 
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3.  Equations  of  motion 

The  equations  are  centered  at 
point  j,k;  see  figure  for  notation. 

I  =  j,k  -  1 

II  =  j  +  l,k 
HI  =  j,k  +  1 
IV  =  j  -  1,  k. 


(a) 


.  n+l/2  _  •  n-  l/2 
Xj,k  =  Xj,k 


At 


n 


2<]>n 


+ 


^  n  n  \ 

-  yjv) 


-B-3- 


UCRL-7322 


n+1 

n  1 

.n+1/2 

x.  . 
J>k 

=  x.  ,  + 
Jsk 

X.  . 

J>k 

n+1 

n  1 

.n+1/2 

yj,k 

=  yj,k  + 

yj,k 

(e)  Calculate  Vq+1  from  Eq.  (2)  (conservation  of  mass). 

(f)  Calculate  q^1^  from  the  equation  of  state  section. 


4.  Equation  of  state 
(A)  Strains 

A®  =  (Aa}®  +  (Ab}®  ’ 

A n+l/2  _  ^  4. 

A®  -  2  r®  +  A® 


u)  (^® 1/2 = (C1/2 = ^ |Az  ■ i4>,y3  ■ yi>  ■ (yz  ■ y4)(i3X_  ii)]n+1/2- 

<*>  Vn1/2  =  fe)r/2  '-T^vz  ih  -  ^  -  X1>  -  (x2  -  ^  -  *>,]n+1/2- 

y  CP  \  /  w 


(Hi)  <Wn+l/2  ■  (I) 


,\n+l/2 


(iv)  (J  ,n+i/2  Jai^r172 
xy  ® 


® 

ax  '  FyyQ 


l® 

T7  -  (e  +  6  ) 

V  xx  yy 


® 


1 


2A 


n+ 

'CD 


{[(y2  "  y4)(y3  "  yl*  "  (y2  “  y4)(y3  "  yl3 
-  [(x2  -  x4)(x3  “  xj)  -  ^x2  x4)(x3  "  x]))| 


(v)  [Ae 


p+l/2  _  \  n+  l/2  Afcn+  l/2 


xx 


Ae 


yy 


'® 

n+1/2 

<D 


XX, 


7® 


yy/ 


n+l/2  Atn+l/2 

® 


Ae 


00, 


n+l/2  =  /.  \n+l/2  Atn+l/2 


Ae 


Xy 


® 

n+  l/2 

® 


ee 


7® 


At 


Xy, 


\n+ 1/2  Atn+ 1/2 

® 


W®  vi  v“+l/z 
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where 


n+l/2  _  1  (Vn+I  ,  vn\ 

v®  -  2 1V(D  +yo  r 


(B)  Stresses 
lo  Elastic 

(i)  u  W1-  (*  \n  +  2l.  Ln+!/2  ifev 
l  ”7®  UU  **  3VV. 


n+l/2 


+  6 


(ii)  (s  'ln+1=fs  \"+2^A,n+l/2-i(^ 

\ yy7®  l ”7®  [  yy  3VV, 


n+ 1/2 


®+  N® 


I  \n,_  f.  n+l/2  1  /AV\n+1^2 

V99)®  Te"  nv)  ® 


(iv)  T 


T  \n  +  +  (S 

i  xy/m  l  xy  \  xy/m 


See  Section  6  for  the  calculation  of  6  ,6  ,  and  6  »  After  the 

n+ 1  XX  Yy  Xy 

stresses  at  time  t  are  calculated,  the  yield  calculations  are  made  (see 

Section  7) . 

2..  Hydrodynamic 

p®  1=a(v®+1) +b(v®+1)-  e®+1  • 

nd  1 

Eq  is  calculated  from  the  total  internal  energy  equation,  Section  5. 

3o  Artificial  viscosity 

n+1/2  co p  A  M2 

(1)  quadratic  q:  q^  -  —~^T/Z -  ‘  l  V  J 


C  2  -  4 
C0  "  4‘ 


Calculate  only  for  V/V  <  0. 


/  „  0,/.  n+l/2 

n+l/2  /^CLP  V 

(11)  linear  q:  q^  =(-  '  -  ^T/T -  V 


=  constant.  Calculate  only  for  V/V  <  0. 


a  =  sound  speed. 
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r 


(qxx)( 


n+l/2  _  [aCA  PW/(A 


<D 


0  A  .n+1/2. 


V 


n+ 1  /2 


(iii)  anisotropic  q:  < 


(V) 


n+l/2  _ 


© 


aCAp°/(A"+1/2) 


rn+  1  /2 


xx 


yy 


<D  ’ 


C  .  =  constant. 
A 


If  these  terms  are  used,  they  are  added  to  the  appropriate  stresses, 

i.  e.  ,  (s  +  a  )  and  (s  +  q  ).  It  is  useful  to  have  the  q  formulated  in 
'  xx  Hxx  yy  ^yy 

this  way  for  problems  where  a  shock  is  travelling  perpendicular  to  a  free 
surface.  However,  for  the  average  problem  the  quadratic  q  gives  very 
satisfactory  results  and  it  is  this  q  that  is  expressed  in  the  equations  written 
here. 

Total  stresses 


4. 

(i) 

(ii) 

(iii) 


- 

n+l 

n+l  r 

2 

xx  _ 

CD  ' 

Sxx 

©  ■  L 

r 

n+l 

1 

n+1  r 

E 

yy-J 

©  = 

s 

yy  J 

©  "  ' 

L  1 

n+l 

r  ' 

n+l  r 

1  b 

CD 

CD 

1  _ 

©  ' 

CD 

CD 

©  ■  L 

n+l/2' 


©* 


|p"+1 

|pn+.  +  y+1/^. 


5.  Energy  equations 
(i)  total  internal: 


TTn+1  -  < 

e(D  ■  < 


|e"  -  [ 


A(Vn+1)  +  Pn  ,  n+l/2" 

2 -  +  q  J 


(yn+1  _  y11)  +  AZn+1//2 


(ii)  anisotropic: 


1  +  E(VT‘I  (V”U  .  V") 

n+1  _  n  n+l/2 

Z®  "  Z  ®+  V® 


<D 


s  e  +  s  e 
xx  xx  yy  yy 


+  s  A  +  T  € 

00  00  xy  xy 


n+  l/2 
© 


AZn+!/2  =  zn+l  _  zn.  sn+l/2  =  ^(gn+l  +  gn}  > 


etc. 

n+1/2 


6.  Correction  for  rotation  of  stresses  during  a  time  step  At 
If  a  mass  element  has  rotated  in  the  x-y  plane  by  angle  go 
during  the  time  interval  Atn  =  t  -  t  ,  the  stresses  must  be  recalculated  so 
that  they  will  be  referred  to  the  x-y  coordinate  system  in  their  new  position. 
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'j’he  following  transformation  equations  can  be  obtained  from  Ref.  2* 


p.  14: 


s 1  “ 

n 

:  g 

2 

cos 

(O 

+ 

n 

s 

.  2 
sm 

CO 

+  2T 

n 

sin 

(O 

cos 

(0 

XX 

XX 

yy 

xy 

s 1  = 

n 
:  s 

.  2 
sm 

CO 

+ 

n 

s 

2 

cos 

CO 

-  2T 

n 

sin 

(O 

cos 

(0 

yy 

XX 

yy 

xy 

n 

r 

2 

2  l 

f  n 

n  1 

OS  (0 

s 

T 1  = 

__  m1A 

cos 

(0 

-  sm  co 

- 

s 

-  j 

3 

c 

xy 

xy 

L 

j 

L  xx 

yyj 

The  angle  co  is  given  by: 

Vx  (ii  -  yj)  =  (||  "  ffH  =  2 


n+l/2 


sin  oo  = 


dy_  _  9x\ 
9x  9y  J 


Equations  (i)  can  be  rewritten  as: 


n  .  n  n  n 

s  +  s  s  -  s 

q.  -  XX  yjr  _jcx _ YJ. 

Sxx  "  2  2 


cos  2(0  +  T  sin  2co 
XY 


n  +  n  sn  _  sn 

-2L- - YIL  _  .  . - — cos  2 co  -  Tn  sin  2co 

2  2  xy 


T1  =  T  cos  2co  - 
xy  xy 


n  n 

s  -  s 
xx  y] 

2 


-  sin  2co 


In  the  incremental  stres  s  -  strain  relations 
n+l  „  .  .  f  1  AV|n+l/2  .. 


n+l 

n 

=  s 

+  2p 

€ 

XX 

XX 

XX 

n 

S  ; 

n 

*  s  , 

yy 

and 

Tn 

XX 

xy 

xx  "  3  V 


,  etc. 


s’  ,  and  T1  .  In 


xx  yy 


additive  term,  S5  to  the  stresses  such  that 

r  .  ...Tn+l/2 

n+l  n  ,  0  1  AV 

S  S  »  2u  €  __  “  0  -rr 


n+  1 

n 

=  s 

+  2\± 

_  ]_ 

6xx  3 

’  /.“ 

XX 

xx 

n  _ 

n 

/  XX 

s  1 

s 

= 

XX 

XX 

XX 

V  ■ 

n 

s'  - 
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7.  Yield  calculations 

(a)  Principal  stresses  (Ref.  5,  p.  94) 


n+ 1  n+ 1 

*3  “  see  ’ 


(In  the  plane  x-y  and  cylindrical  coordinate  systems  used  here  s^g  is  already 
a  principal  stress.  ) 

(b)  Von  Mises  yield  condition 

p  o  lml  afp 

U>  2J"+I  -  [(sf!)2  +  (?f’)2  +  «!)2] 

=  MC/]  K^1)2 

(ii)  2Jn+1  -  2/3  (Y°)2  =  Kn+1  . 

If  Kn+ *  >0  then  multiply  each  of  the  stresses  sn  h  s**  "  ,  and  by 

0  /  J +1 —  +1  xx  yy  at)  xy 

yih  YU//2Jn+1.  KKnTi<  0  use  the  stresses  as  they  are  for  the  next 

time  step. 

8.  Boundary  conditions 

(a)  Fixed  boundary  on  the  x  axis 

Phantom  zones  are  created  by  a  mirror  reflection  across  the 
boundary  as  shown  in  the  figure. 
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The  point  j,  k  can  now  be  accelerated  with  the  equation  of  motion  for 
a  general  point  [Eqs.  (3)]  subject  to  the  following  conditions: 


The  above  procedure  gives  the  desired  acceleration  along  the  boundary, 
but  is  has  the  undesirable  feature  of  not  allowing  for  the  situation  when  the 


point  j,k  is  on  a  free  surface  since  the  point  has  the  extra  mass  of  the 


reflected  zones  associated  with  it.  It  is  more  convenient  to  have  the  correct 

mass  associated  with  each  point,  determined  once  and  for  all  when  the 

problem  grid  is  generated,  and  use  different,  acceleration  routines  for  the 

case  when  the  boundary  is  fixed.  Therefore,  referring  to  the  figure  above, 

we  will  calculate  <j) .  ,  as: 

J  > K 


(ii)  4> 


n 


j  >  k 


/0  An\  /0  n\ 

(£— A--)  +  (P— A_ 
V  vn  /(D  V  V 


n 


The  acceleration  equation  for  point  j,k  that  gives  the  same  results  as 
the  equations  of  motion  for  a  general  point  with  ^conditions  (i)  becomes: 
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(b)  Fixed  boundary  on  the  y  axis 


© 

© 


Similar  to  paragraph  (a)  above,  the  effect  of  a  reflection  about  the  y 
axis,  subject  to  the  conditions  (i) ,  is  obtained  using  the  following  equations 
for  the  acceleration  of  point  j,k: 

f(p°  An/Vn)Q+  (p°  An/Vn)0] 


(ii)  < 


<t>n  V 
\J>k 


pn ,  =  h 

J  »k  2 


'^yy  ‘  ^ee)  (An/M)]  ®+  [C^Vy  "  ^e)  (a  /m)]  <DJ  ' 

(iii)  [dy/dt]  j>k  -  l/2^  k  teyy)®-^)®]  [XII  -  xlll]  "  (Txy)©^II  "  ym] 
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(c)  Corner  zone  on  the  x  axis 


(d)  Corner  zone  on  the  y  axis 

♦ 


Free  surface 
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(e)  Free  surfaces 

For  a  free  surface  at  j,k  in  the  figure  below,  all  quantities  associated 
with  the  phantom  zones  (1)  and  (4)  are  taken  as  zero.  The  equations  of 
motion  for  a  general  point  can  then  be  used,  except  that  aj  }  p.  and  are 

calculated  as  shown  below. 


x 
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(1), 


For  a  corner  free  surface,  the  phantom  zones  in  the  figures  are  zones 
(2),  and  (4).  Free  surface 


Free  surface 


n 

a..  , 
J>k 


p l k 


\  xy  An/M  ® 


2**  -  2 A A 

yy  09. 


An/M 


9-  Sliding  interfaces 

When  two  materials  slide  on  one  another,  a  decoupling  of  grid  points  on 
the  interface  must  be  provided,  otherwise  large  grid  deformations  will  result. 
If  all  of  the  forces  on  the  interface  are  taken  into  account,  the  equations  of 
motion  become  excessively  complicated.  For  a  large  class  of  problems  a 
simple  decoupling  of  the  grid  points  gives  very  satisfactory  results.  The 
method  adapted  here  considers  one  surface  to  be  a  fixed  boundary  during  a 
given  time  step  At.  The  equations  of  motion  for  the  sliding  material  are  the 
same  as  those  given  for  motion  along  a  fixed  boundary.  The  fixed  boundary 
is  then  advanced  in  time  using  the  force  field  of  the  sliding  material  next  to 
it.  The  new  position  of  the  fixed  boundary  provides  a  new  boundary  for  the 
sliding  material.  The  important  point  in  this  type  of  calculation  is  that  the 
parameters  on  the  interface  of  the  sliding  material  be  associated  only  with 
the  sliding  material  and  that  the  material  providing  the  fixed  boundary  be 
treated  as  though  the  boundary  were  an  exterior  surface  with  pressure  forces 
acting  on  it. 
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The  point  is  that  no  forces  should  be  defined  at  the  sliding  interface  by 
an  averaging  process  that  uses  information  from  either  side  of  the  interface. 


Referring  to  the  figure  above: 

(1)  Point  f  is  advanced  along  (k  )n  as  though  (k  )n  were  a  fixed 

s  s 

surface.  The  mass  of  point  f  is  associated  with  the  mass  of  the 

material  below  k  . 

s 

(Z)  Points  a,  b,  etc.  ,  on  (kg)  are  advanced  from  time  (n)  to  time 
(n  +  1). 

These  points  are  associated  with  the  mass  of  the  grid  above  kg. 

The  line  (kg)n  is  considered  a  grid  boundary  for  the  material  above 
and  accelerated  by  forces  from  the  grid  below. 

(3)  Points  on  (kg  -  l)n  are  advanced  in  the  usual  manner. 

(4)  The  point  fn  +  1  if  found  from  the  line  through  Pn+  and  point  *f. 
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Calculations  for  a  sliding  interface 


k  =  slide  line 


(a)  To  calculate  the  volumes  of  k  -  l/2  mass  points 

1.  Given  the  point  P  on  k  -  1  and  the  slope,  m  ,  of  a  line  through  P; 
we  want  to  find  the  points  a  and  b  on  k  that  lie  to  either  side  of  this  line 
(see  figure  )*  Once  these  points  are  found,  the  intersection  f  of  a  line 
through  points  a,  b,  and  the  line  through  P  can  be  determined. 

For  each  point  j,  k  calculate0 

(tan  fl).  k  =  (mp  -  (1  +  mp  •  m._k> 


where: 

m.  .  =  (y.  ,  -  y  )/  (x.  .  -  x  ) 

j,k  '7j),k  'p/  J,k  p 

(x  ,  y^)  and  rn^  are  given. 


(i)  H 
(ii)  If 

(iii)  If 


m 


m. 


m 


>  10^,  (tan  0).  ,  =  l/m. 

J.k  /  j, 

>  10^,  (tan  0) .  ^  =  -  l/m 


>  10  and 


m.  . 

J.k 


>  10  ,  (tan  0)  =  0. 

J»k 
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Test  the  consecutive  values  of  (tan  0).  ,  until  a  change  in  sign  is  found. 

J »  k 

The  points  (x.  .  ,  y.  )  and  (x.  .  ,  ,  y.  ,  ,  )  where  this  occurs  will  be  the 
j,  k  j,  k  j+1,k  j+i,k  .  .11 

points  (xa>  Ya)  and  (x^j  y^)*  (^his  method  fails  if  |0|>  90°.  However,  in 
practice  information  on  neighboring  points  is  carried  from  cycle  to  cycle  and 
these  points  are  tested  first.  If  these  fail,  points  adjoining  on  either  side  are 
alternately  tested.  For  the  original  search,  the  whole  k  line  is  tested  for  a 
change  in  sign  of  arctan  0.) 

2.  To  find  the  coordinates  of  point  f  on  line  ab: 

x  =  (y  -  y  +  m  -  m  ,  x  )/ (m  -  m  ) 
f  wa  'p  p  ab  a'/'  p  ab' 

yf  =(mp  [ya  •  mab  <xa  ‘  Xp)]  ‘  mab  yp)/(mab  ‘  mp> 

mab  =  (ya  -  yb>/<Xa  ‘  V' 

3.  Repeat  steps  1  and  2  for  point  G  on  k  -  1 

4o  To  calculate  the  volume  enclosed  by  Ps  f9  b,  and  G  (see  figure). 
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Aj  =  area  of  Al,  etc. 


2Aj  = 


x,  -  x 

b  p 


xr  -  X 

f  p 


2  A. 


X  -  X 

r  r 


x,  -  x 

b  p 


2A3 


XG  ”  xp 


xr(  “  X 

f'  p 


V 


p+l/2  3 


'S) 


yb  -  yp 

Yf'Yp 


J 

Yf~  Yr 


yb  ~  Y p 


yG  -  Yr 


Yf  ■  yp 


P+1/2 


k 

0 


k 

0 


k 

0 


=  {xb  ~  xp)(yf  "  V  ‘  (yb  ‘  yp)(xf  "  v 


(xf-  -  xP)(yb  -  V  -  (yf  -  yP)(xb  -  v 


=  <*G  -  x P)(yf-  •  y p)  "  (yG“  V(xf-  "  V 


(yP  +  yf  +  yb)Ai  +  (yP  +  yb  +  yf-}  A2  +  {yP  +  yf-  +  yG)j A: 


5.  To  calculate  the  coordinates  of  the  center  of  a  zone 
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Xj+l/2,k-l/2  "  4  (xj,k-l  +  Xf,  k  +  Xf',k  +  Xj+l,k-l) 

yj+l/2,k-l/2  =  4  (yjsk-l  +  yf  ,k  +  yf'3k  +  yj+l,k-l) 

(Note;  This  is  not  in  general  the  center  of  a  zone,  but  is  suitable  for 
the  use  that  is  to  be  made  of  it.) 

(b)  To  advance  in  time  point  f  on  the  slide  line  k 


cos  ab  = 
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I.  Given  the  point  2  in  the  above  figure  and  the  slope  m2  of  a  line 
through  point  2  and  1  to  the  face  ab.  We  want  to  find  the  points  r  and  s 
on  k  -  l/2  that  lie  to  both  sides  of  this  line. 

For  each  point  (j  +  l/2),  (k  -  l/ 2)  calculate: 

(tan  ^>j+ 1/2, k- 1/2  =  (m2  -  mj+l/2,k-l/2>/(1  +  m2  -mj+l/2,k-lA> 

mj+l/2,  k-l/2  =  (yj+  1/2,  k-  1/2  '  y2)A>l/2,k-l/2  -  V2> 

m2  =  -  (xa  -  xb»Aa  -  V 

(.)  I£  mj+l/2,k-l/2  >  '°4’  ,tan  9)j+l/2,k-l/2  =  'l/m2 


(tan  0) . 


j+l/2, k- l/2 


=  -  l/r 


Vl/2, 


^+1/2, k- 1/2 


(tan  0)j+1/2jk_1/2  “  °° 


Test  consecutive  values  of  (tan  0)  j+ 1/2,  k- l/2  until  a  change  in  sign  is 
found.  The  points  (^j  _  2/2,  k- l/2-  l/2,  k- l/2>  and  <xj+  l/2,  k- l/2  * 

yj+l/2,k-l/2>  win  be  the  p°ints  (xr,  yr)  and  (xs,  ys)  (see  details  in  Part  I). 

2.  Calculate  the  point  of  intersection  (xds  yd)  of  the  line  through 
Point  2  and  the  line  through  points  (xr,  yr)  and  (xg)  yg)„ 

Xd  =  (mrsXr  “  m2X2  +  y2  "  yvV^mrs  ~  m2> 

Yd  =  {-2[yr  -  mrs  (xr  -  x2)]  -  mrsy2J/(m2  -  m^) 

mrs  =  <yr  -  ys>/<xr  -  xs> 

If  mrs  <  10"4’  and  m2  >  1()4’  xd  =  x2  and  Yd  =  y  • 

3o  Calculate  the  pressure  at  point  d 
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4.  Repeat  1  through  3  to  get  Pg 

5.  Advance  point  (j,k) 


.  n+ 1/2  .n-l/2  Atn 

x.  1  =  x.  ,  +  - 

j,k  J,k  2<|)n 


Dn  /  n 

pi  y^+ 


1  lyj+3  ,k  yj  ,k+l/  ~2l'j,k+l  7j-],k 


,  Dn  n  n 

+  P9  Y;  V4-1  “  Y^_ 


j>k  +  P^  y 


dl'j-l.k  *j,k  e  r  j ,  k  'j+l,k 


'  n  n 

Yi  v  -  y^+ 


.n+1/2  _  .n-1/2  Atr 


yj,k  =  yj,k 


+  P*  xn 


llj+l.k  j , k+ 1 /  2\j,k+l  j-l,k 


j,k  p.k 


d\j-l,k  j,k/  e\  j,k  j+l,k 


+  pn  xn  ,  -  xn, 


(d)  The  rest  of  the  grid  is  advanced  with  the  normal  equations.  The 
slope  m  of  Section  I  is  found  from  (xf  ,yf)  of  Section  II  and  the  advanced 
position  of  point  P.  We  are  now  ready  to  repeat  steps  (a)  through  (b).  The 
original  slope  m  is  obtained  by  bisecting  the  angle  made  by  point  P  and  its 
neighboring  points  on  line  (kg  -  1). 

(e)  Discussion: 

From  the  above  procedure  it  is  seen  that  the  sliding  material  is  made 
to  follow  the  motion  of  the  slide  line  boundary.  An  error  in  time  and  position 
is  introduced  since,  even  though  mass  is  conserved,  one-half  the  mass  of 
the  sliding  material  is  not  used  in  the  calculation  of  the  acceleration  in  the 
direction  of  the  boundary  motion.  This  error  is  reduced  when  the  sliding 
material  has  a  smaller  density  than  the  material  associated  with  the  boundary 
or  when  small  grid  zones  are  used.  For  reasonable  zoning,  the  error  only 
shows  up  after  large  displacements  have  taken  place.  The  error  can  be 
effectively  remedied  by  increasing  the  mass  in  the  boundary  material  to 
compensate  for  the  one-half  mass  in  the  sliding  material  that  has  been 
neglected  in  the  acceleration  equations. 
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10.  High- explosive  burn  options 

The  chemical  energy  to  be  released  through  the  hydrodynamic  equation 
of  state  is  stored  in  each  high- explosive  zone  as  an  initial  energy  E  .  The 
time  for  the  detonation  front  to  reach  a  specific  zone  is  calculated  in  advance 
from  the  known  detonation  velocity,  D,  and  the  distance  from  the  point  of 
detonation  to  the  center  of  the  zone.  This  quantity  is  also  stored  with  each 
zone.  A  Burn  fraction,  Fs  is  calculated  so  as  to  spread  the  burn  front  over 
several  zones  analogous  to  the  artificial  viscosity  "q"  that  spreads  a  shock 
over  several  zones.  The  burn  fraction  is  integrated  with  the  energy  equation 
and  the  explicit  form  shown  in  Section  5  for  calculating  the  energy  cannot  be 
used.  Instead,  the  energy  must  be  calculated  by  an  iteration  procedure. 

The  procedure  given  below  can  also  be  used  to  replace  the  calculation  shown 
in  Section  5  when  a  more  complicated  equations  of  state  is  required. 

For  one-dimensional  calculations  the  burn  fraction  can  be  defined  as 

F  -  (l-V)/(V  )  (Ref.  11)  and  the  burn  calculation  is  started  by  setting 

F  =  1  in  the  zone  that  corresponds  to  the  point  of  detonation-  The  burn 

calculation  will  proceed  to  around  three  or  four  times  the  number  of  zones 

that  the  artificial  viscosity  nqn  is  spread  over  before  the  detonation  front  is 

correctly  established.  For  one-dimensional  calculations  this  amounts  to 

about  16  zones.  In  two-dimensional  calculations  where  there  is  a  limit  to 

the  number  of  zones  for  a  practical  problem,  it  is  usually  necessary  to  have 

the  correct  detonation  velocity  established  in  a  fewer  number  of  zones.  A 

convenient  way  to  do  this  is  to  start  the  burn  calculation  at  the  time  the 

detonation  would  reach  a  given  zone  as  described  above.  To  allow  for  the 

possibility  of  an  overdriven  detonation  that  may  arise  during  the  calculation 

and  result  in  a  higher  than  normal  detonation  velocity,  the  burn  fraction 

F  =  (1-V)/(1-V^T)  can  be  used  in  addition  to  the  burn  fraction  that  is  based 
J 

on  the  known  detonation  velocity.  The  larger  of  the  two  is  then  selected  for 


the  calculation. 
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(tn+1  -  tb)/AL 


=  (I-Vn+J)/(]-VCJ) 

«n+  1  .  r  T-n+  i 

F  =  maximum  oi  Jt*  ^ 

T  *  -_n+  1  .  ,  ,  ^n+  1  , 

If  F  >  1  set  F  =1 


and  F„ 


(i)  Burn  fraction  =  F^+ ^  -  (tn  -  t^)/AL  AL  -  r  Ax/ D 

,  n+  1  ^  ^  T^n+1  A  t  =  actual  time 

For  t  <  t.  :  ±*  ,  =  U 

t,  =  time  for  a  zone  to 

F^+  =  (I  ~Vn+1)/(l-VCJ)  start  burning 

„n+l  .  ,  ^n+I  ^n+1  Ax  =  grid  spacing 

F  =  maximum  of  F,  and  F^ 

D  =  detonation  velocity 

T  *  -_n+ 1  .  ,  ,  ^n+ 1  _  , 

se  r  =  constant  ~  2.5 

j  =  Chapman- Jouguet 
relative  volume 

(ii)  e"+1  =  E”  -  (Pn  +  qn+l/2)  •  (Vn+1  -  Vn). 

(iii)  Pn+I  =  P(En+1,Vn+1)  •  F1' '  1  , 

here  P(E.,V)  is  the  equation  of  state  of  the  burned  explosive. 

(iv)  En+1  =  En+1  -  l/2(Pn+1  -  Pn)  •  (Vn+1  -  Vn). 

(v)  Pn+1  =  P(En+1  ,  Vn+1)  ■  Fn+1  . 

11.  Stability 

The  At  calculation  is  the  same  as  that  given  for  the  one-dimensional 
problem  in  Appendix  A.  The  characteristic  zone  thickness  is  taken  to  be  the 
zone  area  divided  by  the  longest  diagonal. 

12.  Plane  geometry 

For  plane  geometry  in  x-y  space,  the  mass  calculation  corresponding 
to  Eq.  (la)  at  the  beginning  of  the  appendix  becomes: 


M®=e°/VAa  +  Ab°>- 

The  conservation  of  mass  Eq.  (2)  becomes: 

v®=  (p°/m®)-  (ms)- 

In  the  equations  of  motion  [Eq.  (3)]  the  terms  a  and  (3  are  set  to  zero. 

The  quantity  cj>.  .  is  seen  to  be  a  constant  and  is  calculated  only  once  for  each 
J  *  k 

point  j , k.  The  logic  for  the  value  of  cj>.  ,  at  grid  boundaries  is  the  same  as 

J  s  k 

for  the  cylindrical  case  where  [pb(Ab  +  A^j/V  ]q^  =  »  etc. 

The  term  in  the  equation  of  state  section  is  set  to  zero. 

0  u 
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